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Abstract 


The    Beta-Gamma   transformation    is   described   and    is   used   to 
simple   first-order   autor egress i ve    Beta-Gamma    process,    BGAR(1).      M 
likelihood    estimation    is    discussed    for    this    model,    as    well    as    men; 
estimators.      The    first-order   structure    is   extended   to   include   mi 
processes    and    mixed    first-order    aut or egr essi  ve  ,    pth-ordor   movi: 
processes.      It    is   shown   that    these  Gamma   processes   are   time-revers i I 
therefore,    too   narrow   for   general    physical   modelling.      A   dual    process    to   the 
BGAR(l)    process,    DBGAR(l),    is    introduced,    as    well    as    an    it<  i 
which    combines    the    Beta-Gamma    process   and   the   GAR(1)    process   of   Gaver 
Lewis    (1980).    Some   properties   of    these   extended  autoregressi ve   processes   are 
derived.    Several    highly    nonlinear    extensions    of    these    processes    which 
produce   negative   correlation   are    given.       Use    of    the    processes    to   mode: 
sequence   of    times    between   failures   of   a   computer   system   is   described. 

0.       INTRODUCTION 


The    Gamma    distribution    is    used    to   model    a    wide    variety   of   positive 
valued   random    quantities    in   fields   such   as   operations   analysis,    reliability 
theory,    hydrology    and    meteorology.       Thus,    service   time   distributions 
interarrival    times    in   queues   are  often   modelled    as    having    Gamma    distribu- 
tions,   as    are    wind    velocities    (Hugus,    1982;    Brown,    Katz  and   Murphy,    19 
measured   at    successive   discrete    time    points    or    river    flows    at    success 
instants    of    time    (Lawrance    and    Kottegoda,    1977).      In   all    these   cases,    the 
measurements    are    taken    serially     in    time    and    are    apt    to    t  .    illy 

dependent.      Thus,    development   of    time  series   with   Gamma   distribu 
distributions   and   various   correlation   structures    is   of   great    import 

Gaver    and    Lewis     (1980)    showed    that    the    usual    linear    first-or 

autor  egr  ess  i  ve    equation,    X      =    pX  +E    ,     would    yield    Gamma    mar  - 

&  M  n  n- 1        n 

distributions    if    the    i.i.d.    sequence    {E    }    was   chosen   with   suitable  marginal 

n 

distribution.   This  Gamma  innovation  distribution  has  a  positive  probability 


of  being  zero,  so  that  the  process  (GAR(1))  generates  sample  paths  which 

exhibit  'runs-down'  (as  seen  in  river  flow  data),  but  which  are  "defective". 

The  defect  lies  in  the  fact  that  when  E   =  0,  X   and  X   ,  are  proportional 

n       n      n-1 

and  p  can  be  estimated  exactly  in  long  enough  time  series.  Moreover,  the 
probability  of  the  defect  is  higher  for  k  small,  which  is  precisely  where 
the  model  is  needed,  since  the  usual  techniques  of  transforming  to  normality 
are  then  questionable  and  probably  undesirable. 

Bernier  (1970)  introduced  the  GAR(1)  model  in  a  hydrological  context 
and  McKenzie  (1982)  introduced  a  multiplicative  Gamma  process  called  PAR( 1 ) 
-  product  autoregression  of  order  one. 

In  Lewis  (1983)  and  Hugus  (1982),  a  simple  linear,  random  coefficient 
model  called  BGAR(1)  was  introduced.  It  is  based  on  the  Beta-Gamma 
transformation  described  in  Section  1. 

The  purpose  of  this  paper  is  to  develop  the  properties  of  this  BGAR(1) 
model  and  to  extend  the  idea  to  moving  average  and  mixed  autoregressi ve 
structures.  In  particular,  it  is  shown  that  these  processes,  like  the 
Gaussian  ARMA(p,q)  models,  are  time  reversible  and  therefore  are  very 
particular . 

Several  schemes  for  broadening  the  structure  of  Gamma  time  series  are 
given.  In  particular,  a  technique  of  iteration  produces  a  Gamma  autoregres- 
sive  process  with  two  structional  parameters  that  can  model,  for  given 
marginal  distribution  and  serial  correlation,  different  kinds  of  sample  path 
behavior.  Some  nonlinear  schemes  that  produce  negative  serial  correlation 
are  also  introduced. 

It  is  important  to  note  the  multiplicity  of  Gamma  processes  which  can 
be  derived  with  given  first-  and  second-order  structure.  Conseqently,  in 
the  absence  of  a  'natural'  structure  such  as  exists  for  Gaussian  processes, 
our  aim  has  been  to  produce  simple  structures,  i.e.,  linear,  additive, 
random  coefficient  processes. 

Finally,  a  series  of  times  between  failures  of  digital  computers 
(Lewis,  1964)  is  analyzed  and  fitted  with  the  model.  The  data  is  serially 
correlated  with  a  marginal  distribution  which  is  more  skewed  than  an 
exponential  distribution.  Although  this  data  is  known  to  be  generated  by  a 
branching  Poisson  process  (cluster  process),  the  Gamma  model  is  much  simpler 


and  much  more  tractable  than  the  cluster  process,  and  provides  an  adequate 
representation  for  most  purposes. 

1 ,   PRELIMINARIES 

In  what  follows  we  will  use  B(m,n),  or  simply  B  when  the  parameteriza- 
tion is  clear,  to  stand  for  a  Beta  random  variable  with  parameters  m  >  0  and 
n  >  0,  denoted  by  Beta(m.n).  The  probability  density  function  for  a 
Beta(m.n)  random  variable  is 

fQ(x;m,n)  =  _ ,  ,  _ ,  .  x    (1-x)    ,    0  <  x  S    1  ;  m  >  0  ;  n>  0,      (1.1) 

b  i  { m  )  I  {  n  ) 

where  r(*)  is  the  complete  Gamma  function. 

We  will  denote  by  (B(m,n),  B '  (m' , n '),♦•• }  an  i.i.d.  sequence  of  vector 
random  variables  whose  components  are  independent  Beta  random  variables. 

Let  G(k,6)  stand  for  a  Gamma  random  variable  with  shape  parameter 
k  >  0,  and  rate  parameter  8  >  0,  denoted  by  Gamma(k,B).  The  probability 
density  function  for  a  Gamma(k,8)  random  variable  is 

k  k-1  -Bx 

fG(x;k,6)  =  g  ^(|<)e ,       x  >  0;  8  >  0;  k  >  0.  (1.2) 

We  will  denote  by  {G  (k,8)},  an  i.i.d.  sequence  of  Gamma  variates. 

n 

A    Gamma(k,6)    random   variable   has   moments 
E(G)    =   k/8    =   u;      Var(G)    =   k/82;        C(G)    =   s.d(G/y)    =   k~     2,  (1.3) 

where   C(G)    is   the   coefficient   of   variation,    and   Laplace-Stieltjes    transform 

The  Gamma  variable  is  sometimes  parameterized  in  terms  of  the  parameter 
u  =  E(G).  This  is  useful  in  statistical  work,  since  the  mean  is  a  multi- 
plicative   parameter    and   G   can   be   written   as   a   unit-mean   Gamma   variate,    G*, 


times  u,  i.e.,  G  =  pG*.  However,  in  what  follows,  we  will  use  the  fact  that 
two  independent  Gamma  variates  with  the  same  B-parameters ,  but  possibly 
different   shape   parameters,    add   to   give   another   Gamma   variate 

G"(k+k\B)    =   G(k,B)    +    G'(k',B).  (1.5) 

The  result  is  not  true  if  the  Gamma  variates  have  the  same  mean  but 
different   shape   parameters. 

The   Gamma   family    of    random    variables    include    the    Exponent i al ( k  =  1  )  , 
Erlang(k    integer)    and   Chi~Square(k=r/2 ,    r=1,2,*«*;    6=2)    random   variables. 

Gamma    and    Beta    variates    are    intimately    related    and    two    of    their 
properties   will    be   used   throughout   this   paper. 

(i)   A  Beta(m.n)  variate  may  be  generated  as 

D ,    v      G'  (m,6) / .  m 

B(m'n)  =  G'(m,6)  ♦  G"(n,6)'  (K6) 

where  G'(m,6)  and  G"(n,B)  are  independent.  Furthermore,  the  ratio  B(m,n)  is 
independent  of  the  denominator  G'(m,B)  +  G"(n,6)  =  G(m+n,B)  and  this 
property  characterizes  the  Gamma  random  variable  (Johnson  and  Kotz,  1970a). 

(ii)  The  Beta  Gamma  transformation.  Multiplying  a  G(m+n,B)  random 
variable  by  an  independent  B(m,n)  random  variable  gives  a  G(m,B)  random 
variable 

G(m,B)  =  B(m,n)G(m+n,B) .  (1.7) 

Thus  one  can  reduce  the  shape  parameter  of  a  Gamma  random  variable  (multiply 
by  a  Beta  random  variable),  as  well  as  increase  it  (add  an  independent  Gamma 
variate,  as  at  (1.5)).  A  heuristic  argument  for  the  result  (1.7)  is  that  if 
we  wanted  to  perform  the  operation  (1.7)  on  a  computer,  we  could  first 
generate  G'(m,8)  and  G"(m,8)  to  form  the  ratio  (1.6)  to  obtain  the  B(m,n) 
variate.  There  is,  however,  no  need  to  generate  G(m+n,B)  in  (1.7),  we  can 
use  G'(m,B)  +  G"(n,B),  which  is  independent  of  the  Beta  variate. 
Multiplication  then  gives  G(m,8)  =  G'(m,B). 


(iii)      A    formal    proof   of    the    Beta-Gamma    transformation    is    a    special 
case  of    the   following   Lemma,    which   will    be   used    in   the   sequel. 

Lemma.  Let  X(k,B)  be  a  Gamma  random  variable  and  let  B(kp,kp)  be  a 
Beta  random  variable,  which  is  independent  of  X(k,8),  with  p  =  1  -p 
lying   in    (0,1).      Then 


E{e 


-(v+Bu)X      __    (_^]kp    (_J_]kp 

lB  +  v  '         lB  +  v+u  '  =      '         = 


(1.8) 


When  v  =  0,  this  result  proves  the  Beta-Gamma  transformation  in  the 


form 


X(kp,B)  =  B(kp;kp)X(k,B) 


(1.9) 


Proof .   Using  (1.4)  and  conditioning  on  B,  we  get 


„,    -(v+Bu)X1         „ 

E{e  j    =   Er 


i  +  v+Bu 


6*v     B-|^ 


U0  -*V 


k  •  i  rrifel  «Bi), 


where  we  have  used  the  finiteness  of  the  expectation  to  take  the  expectation 

inside  of  the  binomial  expansion. 

Now  since  E(B')  =  B(  k  p  +  £,  k  p  )/B(  k  p,k  p  )  =  r(k  p  +  D  r(k  p  )  r(  k  )/{  r(  k  +  8.)  r(  k  p  )  r(k  p) 
we  have  that 


k  +  H  L.  i.  rk  +  i-i  ^   r(kp  +  i,)r(kp)r(k) 
1  r(kU)r(kp)r(kp) 


r(k  +  g-)  r(kp  +  2.)  r(kp)r(k) 

r(k)r(  a+1  )r(kp)r(kp)r(k  +  H) 


kp+Ji-1 
i 


Thus 


k  +  fc-1 


U       t       )(-   ^)    MB    ,    -      U       t         )(-   ^ 


kp         ,  Nkprni~kp 


u  6+v+u  8+v 

1  -• 

B+v 


and   therefore 

„,    -(v+Bu)X,         (    B    1kpf_6_1kp 

which   was    to   be   proved. 

2.       THE    FIRST-ORDER  AUTOREGRESSIVE   PROCESS,    BGAR( 1 ) 

2.1.   Construction  of  the  Process. 

Using  the  Beta-Gamma  transformation,  we  can  construct  a  very  simple 
first-order  autogressive  process  {X  (k,8,p)}  with  Gamma(k,8)  marginal 
distribution  and  a  single  parameter,  p,  that  describes  the  dependency 
structure  of    the    process.      We   have 


X    (k,B,p)    =    Bn(kp,kp)Xn_1 (k.B.p)    +    B^(kp,kp)G    (k,B)  (0    <    p    <    1 )     (2.1) 

=    Bn(kp,kp)Xn_1 (k,6,p)    +    Yn(kp,B)  n=0,±1,...,  (2.2) 

where    {Y    (kp,B)}    are    i.i.d.    Gamma(kp,B)    variates    independent    of    the 

{B     (kp.kp)}     sequence.        If    X    _     (k,B,p)     has    a    Gamma(k,8)     marginal 

distribution,    then   multiplying    by    B    (kp,kp)    reduces    it    to    a    Gamma(kp.B) 

variate    and    adding    the    innovation   variable   Y    (kp,B)    creates    the   Gamma(k,6) 

n 

variate,  X  (k,B).  The  alternate  form  (2.1)  shows  the  process  as  a 
transformation  of  an  i.i.d.  Gamma(k,B)  sequence,  but  clearly  generation  on  a 
computer  would  be  done  with  (2.2).  In  the  sequel,  we  will  drop  the 
parametric  notation  where  no  confusion  is  possible. 


It  is  clear  that  taking  X  to  be  a  Gamma(k,B)  variate  will  start  the 
process  in  a  stationary  mode.  Also,  the  process  is  Markovian  by 
construction. 


2.2.      Serial   Correlation. 

It    is   easily   established,    using  moments   of   Beta  variables,    that 

p(r)    =    Corr(X    ,X         )    =    p'r',  r   =   0, ±1 , ±2, • • •    .  (2.3) 

n      n-r 

Thus,  in  the  three  parameter  process,  the  parameters  k  and  8  describe  the 
marginal  distribution  of  the  process  and  p  independently  describes  the 
correlation  structure.  Note  that  since  the  process  is  only  defined  for 
0    ^    p    <    1    the   correlations   are  non-negative. 

2.3-      Joint   Laplace-Stieltjes    transform. 

In    the    stationary    process    {X    },    let    L  (u,v)    denote   the   joint 

n  a    ,  a       , 

n      n-1 

Laplace-Stieltjes    transform   of    the   adjacent   variables   X      and   X      ,  .       Then    we 
v  J  n  n-1 

have 


Lv      v         (u,v)    =    E[exp{-X   u-X      .v}] 
X    ,X      ,  n        n-1 

n      n-1 

=   E[exp{-B   X      ,u-Y   u-X      ,v}] 
n   n-1        n        n-1 

=    E(e-uY)E{e-(v+Bu)X},  (2.K 


where,    in   the   last   step,    we   have   dropped   the    indices    n    and    n-1    because    of 

stationarity   and   have   used   the   assumed    independence   of   Y      and   X  to  write 

J  n  n-1 

the   expectation   as    the   product   of   two  expectations. 

Now   the   second   term   is   evaluated    in   the   Lemma  of  Section    1    and  we   have 


,  ,         s         (    B    ikpr    8    ^kpr      8      )kp 

Lx   ,x     ,(u'v)  "  fe]     W     [JW^] 
n     n-1 


,kpr      6      ikp 


P_     #     _LjKPf_±_|Kp 


Since  this  transform  is  symmetric  in  u  and  v,  the  joint  distribution  of 

X   and  X   .  is  symmetric.   Also,  since  the  joint  distribution  of  any  set  of 

n      n-1     J  J 

X  's  can  be  obtained  from  (2.5)  and  the  marginal  Gamma  distribution,  the 

Beta-Gamma  process  is  time-reversible. 

Note,  too,  that  we  have  directly  from  the  defining  equation  (2.2)  that 

the  regression  of  X   on  X   ,  =  x  is  linear: 

n     n-1 

E(xnlxn-rx)  =  pX  +  n-p)k/B  =  px  +  (1-p)u.  (2-6) 

The  time-reversibility  of  the  process  shows  that 


E(x  _i  lxn=y)  =  Py  +  <  T  ~P)  vj-  (2.7) 


2.4.   Convergence  to  a  Gaussian  AR(1)  Process. 

1  /2 
As  k  gets  large,  the  standardized  Gamma(k,1)  variate  X'  =  (X~k)/k 

converges  to  a  standardized  Gaussian  variate.   To  prove  this,  consider  the 

pair  X   and  X   ,  in  the  BGAR(1)  process.   The  joint  characteristic  function 
n      n~1 

for  the  standardized  variables  X'  and  X'  ,  is,  using  (2.5), 

n  n-1 

4>x,    x,       (s,t)    =   E{exp(-isX]^-itX^_1)} 

n'    n-1  .  _  -  - 

ik1  /d(s+t)l.    .   ,-1/23     PE.    .,.,-1/23"     ,,    .,-1/2,     ,u-kp 
=   e  [1+isk  j        [1+itk  ]         {1+ik  (s+t)} 

-1  /2 
Taking   logarithms   and   expanding    in   powers   of    k  gives 

4»x,    x,       (s,t)    =    ik1/2(s+t)    -    kp{isk~1/2-(is)2/(2k)+0(k~3/2)} 

n'    n_1  -1/2  2  -3/2 

-  kp{ itk         -(it)   /(2k)+0(k   3      )} 

-  kp{ ik~1/2(s+t)-(i)2(s+t)2/(2k)+0(k~3/2)} 


and  as  k^00,  this  converges  to 

Yx,  x,   (s,t)  =  -  |{s2+t2+2pst}.  (2.8) 

n'  n-1 

Thus,  since  the  process  is  Markovian,  the  BGAR(1)  process  is  equivalent  to  a 
Gaussian  AR(1)  process  when  k  becomes  large. 

2.5.   Additivity  and  the  GAR(1)  Process. 

Gaver    and    Lewis    (1980)    showed    that    the    usual     linear    stochastic 
difference   equation 


Xn    ■    PXn_,     *    E„.  (2.9) 

will  give  a  process  (GAR(1))  with  Gamma(k)  marginal  distributions  if  E   is 

n 

the   Gamma    innovation   variable    GI    (k,p)    with    Laplace-St i el t j es    transform 

k  n 

{ ( 1 +ps)/ ( 1 +s) }    ,    where  0   S    p    <    1.      This   variable   can   be   simulated   by   methods 

given   by   Lawrance    (1982)    and  McKenzie    (1986). 

Note    that    the    GAR(1)     process     is    a    linear    additive     (constant 

coefficient)    process   and  adding   two   independent    GAR(l)    processes    with    the 

same    8    and    p    values,    say    (X*(k, )}    and    {X**(k„)}    gives   a  new  GAR(1)    process 

n      1  n        2 

with   shape    parameter    k 1 +k      and   dependency   parameter    p.      This   is   not   true   for 
the    BGAR(1)    process   which    is   a  random   coefficient,    linear   additive   process. 

The   process    obtained    by    addition    is    a    process    with    Gamma    marginals    and 

rl 

correlation  structure   p(r)    =    p'     ',    but    it    is   not    even   Markovian.      Additional 

differences    between   the    processes    are    that    while    the    BGAR(1)    process    is 
time-reversible,    the    GAR(1)    process    exhibits    'runs   down'.      In   fact,    it    is 

degenerate    in    the    sense    that    the    innovation    variable   GI       is    zero   with 

k  \ 

probability    p    .      Thus,    we  get   X    /X  =    p   with   probability    p    ,    and    p   can   be 

estimated    exactly    in    long    enough    series    (Gaver    and    Lewis,     1980).       This 

degenerate    behavior    is   not    exhibited   by   the   BGAR(1)    process.      A  method   for 

combining   two  processes    to  obtain  a  broader   process   is   described    in   Section 

5. 


2.6.   The  Conditional  Density  for  X  ,  Given  X   .  -  y. 
From  the  definition  (2.2),  we  have  that 

P{X   s  x  l  X   ,  =  z}  =  P{B  (kp,kp)z  +  Y  (kp.B)  ^  x} 
n     '  n-1  n  n 

=  P{Y  (kp.B)  ^  x  -  B  (kp,kp)z}. 
n  n 

Now,    by   definition,    Y    (kp,B)    is   independent   of   X  and  of   B    (kp.kp).       Thus 

conditioning    on    B       and    differentiating    with    respect    to    x    yields    the 

conditional    density   for   X    ,    given   X      ,    =   z ,    as 
J  n     °  n-1 

fY    ,v        (x,y)    =  11^— x   kkpe"kX  x   >   0,    y    I   0 

Wi  r(kp){r(kP)}2 


L 

kp-1  , .       . k p-1  .  .    kyw  .  . 

w  (1-w)  (x-yw)e         dw,  (2.10) 


where 

I      1  if    x    ^   y, 

L  =  (2.11) 

{      x/y        if    x    <    y, 

and,    for   simplicity,    the   scale   parameter    6   has   been   set    equal    to  one. 

This  density  is  a  continuous  function  of  x  and  absolutely  continuous 
except  where  x  =  y.  its  utility  is  in  obtaining  maximum  likelihood 
estimates  of  the  parameters  8  (or  m),  k  and  p.  This  is  discussed  in  the 
next   subsection. 

2.7.      Moment    and   Maximum  Likelihood   Estimates    in   BGAR(1). 

There  are  "natural"  moment  estimators  for  the  three  parameters  in  the 
BGAR(  1  )  model,  namely  the  mean,  \i,  (or  8),  the  shape  parameter,  k,  and  the 
first-order  serial  correlation  coefficient,  p.  From  (2.3)  and  (1.3) ,  these 
estimates   are,    from   a   sample    (observed   time  series)    of   length  n, 
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n 
p  =  X  =  I   X  /n,  (2.1?) 

i  =  1 

k  =  (X)2/S2,  (2.13) 

where  S2  is  the  sample  variance,  and 

n-1 
p  =  I    (X  -X)(X    -X)/{(n-1)S2}.  (2.14) 

i  =  1 

n 
The    variance    of    p    is    var(X)[1+2    I    (1-(r/n)}p    ]/n    and    nVar(u)     is 

r  =  1 
asymptotically    equal    to    p  2  (  1  +  p  )  /{  (  1  -p)k  } .      More  general    properties   of    the 

estimates   R   and    p,  however,    are  hard   to  derive.       But,    their    distributions 

are    independent    of  the   scale   parameter    p. 

In  Table   2.1,  we   give   the   simulated   standard   deviation   and   bias    of    the 

estimates      K    and  p    for    values    of    k    =    0.25,    0.75,     1.0,    2.0   and    p   =  0.25, 

0.60,    0.90    for    various   values   of   n.      Note  that    the   values    of    n    differ    with 

1  /2 
p,     since    the    "equivalent    sample    size",    n'=n[(1+p)/{(1-p)k}]         ,    is 

different.      Here   n'    is   the   sample  size  which   would   be    needed,    for    given    p, 

to    achieve    the    same    variance   for    p   as    in   the    p=0    (independence)   case.      The 

simulation   study  was   performed  with   the   SUPER-SIMTBED   program    of    Lewis,    et 

al.    (1985). 

Two  other  properties  of  the  process  may  be  useful  in  validating  the 
BGAR( 1 )    model    from   data. 

The   first    is   that    difference   of   successive   values    in    the    time    series, 

D      =    X      .  -X    ,    have    an    2,-Laplace   distribution    (Dewald,    1985).      This   result 
n  n+ 1       n 

comes   from    (2.4)  , 


l  1 


o 
o 


n 


O 

o 


o 

ON 


O 

II 


IfN 
C\J 


O 

o 


CNJ 


O 
O 


i  J*l 


o 

O 
II 
Q. 


o 


in 


LPl 
CM 

O 

II 

Q. 


O 

o 


CM 


o 
o 


in 

r 


«-  *x> 

oo 

f- 

LPl 

OO 

ON 

CO 

ON    C— 

so 

CO 

, 

LfN 

•=r 

( 

oo  «- 

CNJ 

t— 

OO 

in 

< 

oo 

oo  so 

t-- 

4 — 

SO 

CO 

CNJ 

•a- 

in  «- 

ON 

LT\ 

vO 

oo 

vO 

CNJ 

LTl    CNJ 

^r 

CNJ 

^r 

«- 

•=r 

*— 

,_  ,_ 

O 

o 

o 

o 

O 

o 

O    O 

o 

o 

o 

o 

o 

C  1 

+ 

+ 

+ 

+ 

+ 

+ 

+ 

+ 

LTl  t— 
ON  UD 
CO   so 


O    O 

+ 


o  o 

+ 


CNJ    O 
ON  ON 

oo  cnj 


o  o 

+ 


O    O 

+ 


SO  vO 

o  ^r 
LPi  oo 


ON    CNJ 
OO   CNJ 


CO    CNJ 

^r  o 

OO    CNJ 


.-    CO 

O  oo 

OO   i- 


OO  CNJ 
CO    — 

CNJ    ,- 


c--  in 

\D    CNJ 
CNJ    t- 


O    O 

+ 


O    O 

+ 


O    O 

+ 


o  o 

+ 


o  o 

+ 


o  o 

+ 


oo  in 

CNJ    [^ 

CNJ    CO 

t-  in 

ON    t- 

co  ^r 

oo  .=r 

T-     CO 

CNJ    CO 

oo  m 

OO    CNJ 

CNJ      r- 

t—  m 

•=r   cnj 

oo  .— 

C\J     .- 

CNJ    .- 

CNJ    >- 

ON 

o 


o  o 

+ 


o  O 

+ 


o  O 

+ 


o  o 

+ 


o  o 

+ 


O    O 


CO     ON 

-—  in 

CNJ    ,- 


oo  oo 

r-     O 


in  cnj 

«-   CO 

•-  o 


t~-  SO 

o  ^o 

*-  o 


O    =T 
O  SO 

•-  o 


O    O 

+ 


O    O 

+ 


o  o 

+ 


o  o 

+ 


o  o 

+ 


o  o 

+ 


o  o 

+ 


o 

OO 


o  o 

+ 


o 

CNJ 


o  o 

+ 


o 


O    O 

+ 


o 
o 

CNJ 


o  O 

+ 


o 

CNJ 


o  o 

+ 


o 

CO 
CNJ 


OO    y- 

m  o 

CNJ    ,- 


o  o 

+ 


-=T    ON 
ON    CO 

i-   o 


o  o 

+ 


^r  oo 
on  m 
o  o 


o  o 

+ 


CNJ 

CNJ 

^r 

SO 

n 

SO 

t-- 

CO 

CO 

SO 

in 

=T 

oo 

CNJ 

«— 

CNJ 

OO 

^T 

in 

so 

t~- 

o 

o 

CNJ 

oo 

CO 

ON 

oo 

T_ 

o 

-=r 

OO  vD 

o 

CNJ 

t~- 

oo 

CO 

CNJ 

OO 

CNJ 

in 

CO 

o 

SO 

c— 

in 

in 

^r 

oo 

•=r 

1 — 

CNJ 

so 

CNJ 

=T 

"— 

OO 

o 

oo 

o 

CNJ 

o 

CNJ 

o 

CNJ 

o 

CNJ 

O 

o 

o 

o 

o 

o 

o 

o 

o 

o 

o 

o 

o 

o 

o 

o 

n 

+ 

+ 

+ 

+ 

+ 

+ 

+ 

+ 

=r   cnj 

SO   o 

o 

OO 

ON 

C— 

.— 

in 

in  m 

f—    OO 

oo 

CNJ 

CNJ    .- 

^r   c~- 

oo 

i — 

^r 

ON 

CNJ 

f- 

O   SO 

CO  so 

^o 

in 

in  oo 

oo  .— 

CNJ 

t— 

CNJ 

o 

CNJ 

o 

CM    o 

«-  o 

<— 

o 

O    O 

+ 

O    O 

+ 

o 

o 

+ 

o 

O 

+ 

o 

O 

+ 

o 

o 

+ 

o 

o 

+ 

o 

o 

+ 

OO    CM 

oo  t— 
=r   CM 

OO  vO 
CO    -3" 
CNJ    <- 

o 
oo 

CNJ 

o 

m 
o 

CNJ 

CNJ 
OO 

o 

in 
c— 

OO 
O 

ON 

in 

CNJ 

in 
o 

CNJ 

in 

o 

CNJ 

in 

ON 
OO 

o 

o  o 

+ 

o  o 

+ 

o 

o 

+ 

o 

o 

+ 

o 

o 

+ 

o 

o 

+ 

o 

o 

+ 

o 

o 
+ 

OO  C~- 

«~  in 

CNJ    r- 

in  co 

OO  OO 

t-   o 

^ 

o 
o 

ON 

O 

in 
o 

CO 
O 

OO 

o 

o 

in 

OO 

o 

o 

OO 

o 

o 

ON 

CNJ 
O 

O    O 


o 

CNJ 

oo 


CNJ 
CNJ 

CNJ 

in 

oo 

ON 

CNJ 

in 

CNJ 
CO 

CO 

-3- 

in 

CNJ 

o 

OO 

oo 

ON 

o 

OO 

ON 

t— 
o 

oo  og 
oo  t~- 

oo  o 

' 

+ 

o 

o 

+ 

o 

o 

+ 

o 

o 

+ 

o 

o 

+ 

o 

o 

+ 

o 

o 

+ 

o  o 

+ 

in 

OO 

oo 

in 

ON 
OO 

ON 

t>- 

OO 

o 

CNJ 

CNJ 

vO 

ON 

o 

CNJ 

oo 
o 

CNJ 
CNJ 

o 

o 

CNJ 

o 

O 

i-  t-- 
on  in 
^-  o 

o 

o 

+ 

o 

o 

+ 

o 

O 

+ 

o 

o 

+ 

o 

o 

+ 

o 

+ 

o 

O 

+ 

o  o 

+ 

CNJ 
ON 

vO 
CO 

CNJ 

OO 
OO 

in 

CNJ 

in 
o 

CNJ 

co 

o 

ON 

CNJ 

o 

in 
tr— 

CNJ 

in 
o 

in 
o 

co  [■— 
in  ^r 

«-  o 

o 

o 

+ 

o 

+ 

o 

o 

+ 

o 

o 

+ 

o 

o 

+ 

O 

o 

+ 

o 

+ 

O    O 

+ 

CNJ 
CNJ 

oo 

oo 

ON 

o 

oo 

o 
o 

ON 
ON 

o 

in 
o 

CNJ 
ON 

o 

o 

oo 

o 

o 
o 

UN 

o 

oo 
oo 
o 

OO  ON 
f-    C\J 

o  o 

O 

o 

+ 

o 

o 

+ 

o 

o 

+ 

o 

o 

+ 

o 

o 

+ 

o 

o 

+ 

o 

o 

+ 

O    O 

+ 

>, 

rO 

4-> 

■—I 

c 

4J 

••-t 

c 

CO 

* 

D 

^c 

cr 

m 

CU 

. 

s~ 

x: 

J* 

d) 

4-> 

4-> 

«i_i 

0) 

^H 

o 

ct3 

r-H 

c 

L, 

0) 

o 

CO 

a 

•^H 

a. 

4-> 

x: 

CO 

CD 

o 

•r-l 

Q. 

CO 

> 

CO 

CD 

<D 

JC 

TO 

n 

C 

I— 1 

•o 

0) 

L 

£ 

ffl 

4-> 

• 

•o 

c 

C 

Oh 

CO 

O 

CO 

4-1 

CD 

co 

* 

IS] 

1    J* 

■  rA 

x> 

CO 

CD 

•• 

4-> 

£- 

CD 

CO 

O 

rH 

E 

jJ 

Q. 

•^H 

CO 

F 

4-> 

e 

CO 

co 

•^ 

CO 

CD 

4-> 

cn 

T3 

CD 

CD 

C 

x: 

m 

4-> 

j-> 

c 

m 

CO 

cfl 

CD 

Q. 

'*-* 

• 

e 

1 — 

o 

-o 

>> 

• 

E 

c 

4J 

CNJ 

CO 

■^ 

QJ 

4-> 

UJ 

x: 

X. 

c 

_) 

j-> 

CO 

CO 

co 

3 

<x 

Vh 

i- 

cr 

E- 

o 

CD 

4-> 

Q. 

CO 

CD 

O 

CO 

f= 

4-> 

•l-i 

CO 

J3 

C- 

CD 

CO 

x: 

■D 

a 

4-> 

C 

CO 

CD 

■a 

x: 

c 

c 

4-> 

m 

o 

•^H 

C^H      1 

^ 

4-> 

o 

CO 

U-, 

•f-H 

co 

o 

> 

CD 

OJ 

3 

co 

■o 

rH 

CO 

CO 

•rH 

■a 

> 

JD 

L 

CO 

co 

■a 

T3 

3 

<D 

C 

O 

4-> 

CO 

■rH 

CO 

4-5 

t- 

E 

CO 

CO 

•rH 

> 

4-> 

d) 

CO 

Jd 

tn 

CD 

JJ 

o 

Ci-i 

CD 

4h 

x: 

O 

co 

CO 

4-5 

>. 

CD 

to 

■o 

O 

•rH 

D 

o 

CO 

4-> 

t. 

•rH 

CO 

a. 

CO 

CD 

c 

CO 

x: 

o 

F 

4-> 

•l-i 

F 

c 

4J 

cO 

0 

CO 

a 

5- 

^H 

i 

CO 

3 

(0 

Q. 

e 

4-> 

■r-4 

CD 

C 

CO 

CD 

■rH 

C  r- 


-3" 

OO 


CNJ 

m 


12 


oo 

CO 


o 
o 


oo 

oo 


o 
o> 

o 

II 

0_ 


o 
o 


CM 


o 


LCTO 


in 
c\j 


o  on 

ON 

^T 

CO    .- 

in  oo 

o  m 

00 

CM 

m  cm 

on  o 

[>-  vo 

^r 

on 

on  cm 

on  <— 

<M   .- 

U1 

. 

CM     r- 

CM    — 

O  O 

o 

o 

o  o 

o  o 

o  o 

o 

o 

o  o 

o  o 

cl 


1     Q-l 


O 

o 
It 

Q. 


o 
o 


o 
o 


m 


in 

CM 


cl 


.1  , 

CM 
O 

o 
II 

d 


o 
o 


CM 


: 
O 


in 
t-- 


in 

CM 


o  O 

1 

o 

o 
1 

o 

o 
1 

<  ) 

o 
1 

'  > 

-  ■ 
i 

I 

> 

1 

1 

oo  [-- 
o  o 

7> 

Lf- 

o 

o 
o 

o 

00 
CM 
O 

J=T 
<    1 

f  J 

f  J 
o 

(  i 

in 

• 

CM 
O 

o 

C  > 

00 

!  i 

o  O 
1 

1  > 
1 

<  > 

O 

1 

1 
1 

-  ) 

i 

■ 

1 
I 

1  1 
1 

! 

i—  on 

ON    t— 

o  o 

U  i 

o 

on 
O 

f  1 

IJ . 
o 

a 
C  J 

o 

-r 

i 
i   i 

CM 

o 

r 
O 

■ 

m 

• 
o 

1 

i 

on 

o 

.1 

on 

o 

o 

1  1 

o  o 

1 

o 

1 

o 

o 

1 

C  .) 

o 

! 

' 

• 
I 

o 

o 
1 

1 

o 

on  a\ 

«-  o 

o 
o 

O 

LTO 

.J 
o 

O 

U 
f— 

I 

O 

o 

t- 

m 

o 

CM 

o 

o 

i;  i 
o 

T 

CM 
O 

in 

IT 
O 

li  1 

CO 

o 

o 
If 
o 

O 

< 

o  o 
1 

o 

o 

1 

o 

o 
1 

o 

o 

1 

c .) 

o 
1 

I 

1 

• 

1 

"  - 

1 

CM 

CM 

--J 

kO 

CO 

CO 

-i.) 

m 

r 

CM 

on 

i 

o 

r»  i 

CM 

L> 

[ 

71-, 

OO 

in 

o 

TO. 
CM 

CM 

r.  i 
O 

o 

-T 

CM 

O 

CM 
C3> 

o 

oo 

C  3 

O 

in 
o 

in 

r 

o 

OO 

o 
o 

r 

C3 

o 

o 
o 

o 

o 

o 

1 

O 

o 
1 

o 

o 

1 

o 

O 
1 

( 

o 
1 

1 

1  1 
1 

:  ) 

1 
1 

o 

1 

t-~ 

f'l 

CM 

v£3 

O 

- 

o 

o 

c  > 

in 
m 

o 

CO 

o 

a  i 
CM 

o 

0 
CO 

o 

m 

OO 

O 

o 

r  i 

r  j 
O 

T 

t- 
o 

oo 

o 

o 

1 

O 

o 
1 

o 

o 
1 

o 

o 

1 

o 

o 
1 

o 

o 
1 

o 

O 
1 

: 

1 

I 
co 

o 

CM 

=J 

to 
o 

o 

CM 

o 

LO 

o 

in 

o 

to 
m 
o 

a 
a 

o 

li 

on 

o 

o 

CM 
O 

'I 

o 

in 

CM 
O 

CM 
CO 
O 

CM 

CM 

C 

o 

o 
1 

O 

o 
1 

o 

o 
1 

o 

o 

i 

o 

c 

i 

3 

O 

1 

c 

o 

1 

- 
1 

CM 
CM 

=J 

on 
O 

o 
o 

o 

[■-- 
on 

o 

in 

o 
o 

r>  i 

o 

CM 

t   - 

o 

CM 
CM 

-I 
o 

o 

X 

rt  i 
O 

O 

o 
1 

o 

o 
1 

o 

o 

1 

o 

o 

I 

O 

o 
1 

- 

o 

1 

i    ) 

c 
1 

c 

C  3 

i 

o 

o 

o 

o 

O 

o 

o 

r< 

oo 

i© 

o 

-I 

a; 

CM 

i— 

I— 

O.J 

rj 

Oj 

on  cm 

CM   .— 

in 

CO 

o 

in 

CO 

o 

on  co 
on  on 

r~     O 

OO 

t 

00 

o 

on 

on 

CM 

o 

-I 
o 

C  J 

I 
o 

o  o 

I 

o 

o 

1 

o 

o 

1 

o 

o 
1 

( 

o 

1 

c 

o 
1 

o 
1 

< 

o 

1 

^r    in 
on  on 
CM   .- 

c— 
oo 

o 

<o 

vD 

in 

o 

r 
O 

• 
o 

Oj 

o 

o 

^o 

.1 1 

00 

o 

o 

7> 

CM 
O 

o  o 
I 

o 

C  3 

1 

o 

o 
1 

O 

o 
1 

1 

o 
1 

> 

o 

I 

1 

c 
1 

1 

on  on 

00    ^~ 

a 

o 

in 

00 

in 
o 

O 

in 

c 

o 

o 

J 

o- 
m 

o 

• 

CM 

o 

0 

o 

U  '■ 

■j. 

CM 

O 

o  o 

1 

3 

o 
1. 

O 

C  ' 

i 

o 

o 

1 

(  3 

1 

O 

o 
1 

1 

o 

1 

1 

OO    CO 

oo  ,- 

CO 

O 
CM 

o 

30 

OJ 

o 

o 

o 

-i 

(  i 

in 

J 

o 

CO 

o 

3 

o 

o  o 

1 

( 

o 

1 

O 

o 

I 

O 

o 
1 

C  ) 

o 

i 

<- 

1 

<  3 

o 

1 

' 

o 

1 

- 

*J 

L 

C 

ai 

cfl 

»j 

3 

cu 

CT 

i 

11 

L 

£ 

co 

4J 

• 

a 

1 

Q. 

c 

^■1 

o 

^H 

0 

•  ^ 

tu 

CJ 

.: 

rfl 

0 

r   < 

x: 

.^ 

QJ 

CJ 

i. 

TO 

TO 

C 

■:■ 

■^« 

o 

> 

CJ 

c 

•: 

M 

•a 

QJ 

x: 

tj 

4J 

• 

:. 

C 

TO 

-3 

o 

w 

c 

"- 

TO 

» 

I . 

1  1 

1    o 

v. 

TO 

- 

•o 

t. 

0) 

a 

' 

-  ' 

- 

a 

;•• 

ro 

r 

E 

V. 

rfl 

.-. 

05 

- 

TO 

•/j 

a 

CD 

i; 

c 

CO 

OJ 

i- 

.TO 

Q 

c 

n 

^J 

• 

a 

Q 

^ 

E 

TO' 

• 

O 

■a 

•rH 

CM 

E 

c 

CO 

>. 

u: 

0) 

iJ 

.  J 

c 

.1 

■^H 

:. 

t-1 

i 

•X 

m 

TO 

r- 

4_ 

L 

rfl 

O 

OJ 

TO 

JJ> 

TO- 

■J; 

OJ 

- 

^ 

.- 

■  ^ 

S 

0 

£ 

- 

■D 

a 

- 

C 

c 

~ 

] 
c 

i 

: 

*j 

Z' 

O 

TO 

.^j 

-». 

TO 

♦ 

o 

CO 

Q. 

•r-l 

00 

> 

1 

Cl  . 

a 

3 

c 

'■ 

rH 

TO 

TO 

x- 

> 

TO 

L 

.,  -. 

rfl 

to 

TO: 

" 

3 

C 

O 

TO 

CO 

.^ 

0) 

* 

L 

.  ' 

<:. 

TO 

cfl 

> 

E 

1 

.,H 

C 

L 

4 

4-3 

O 

TO 

1 

I  . 

0 

CAJ 

a 

TO 

TO 

>, 

0) 

*  1 

--■ 

c 

-- 

0 

TO 

* 

t. 

■  ^. 

m 

... 

(TO 

.: 

TO 

a 

O 

E 

TO 

.^H 

t" 

a 

- 

TO 

TO 

CO 

C^ 

r— I 

1 

TO 

3 

CO 

C 

E 

*J 

L 

■•H 

: 

TO 

CO 

DD 

a 

7 

on 


00 

in 


o 


rv 
00 


O 
O 


rv 


13 


after  converting  to  characteristic  functions  and  setting  v  =  -u.   We  have 
the  characteristic  function 


>D  (u)  =  Lx      (iu.-iu) 
n        n  n-1 


;  +  iu   8-iu 
„2  ikp 


kp  r  n  -\ k p 


6  +u 


This   is   the   characteristic   function  of   an    St-Laplace   random   variable  with 

2.   =   kp;    the    distribution    is    symmetric    about    zero.       It    goes    to    a    Normal 

random    variable    as    i    ■*    °°,    but    for    I    S    1  ,    the    density    function    is    not 

absolutely    continuous    at    zero.       In    fact,    for    I    ^    1/2,    the    density    is 

infinite    at    zero.       The    fact    that    the    difference    has    median    value    0, 

irrespective   of    the   value   of    p   or    8,    can   be   useful    in   validating   the  model. 

Consider,    now,    ratios    R      =    X      , /X    ;    from    (2.2)    this    is 

n    n+1   n 


R   =  B 
n 


n+1(kp,kp)  ♦  Yn+1(kp,8)/Xn 


Bu 


t  Y   ,(kp,B)  and  X   are  independent  Gamma  variates,  so  that,  if  k>1 , 
n+1  n         v 

E(R  )  =  E{B   .(kp.kp)  +  E{Y   .(kp,8)/X  } 
n       n+1  n+1        n 

=  p  +  (kp/8)/{(k-1)/6}  =  1  +  p/(k-1).  (2.15) 

Higher  moments   can   also   be   obtained. 

These  results  could  be  of  use  in  validating  the  model.  Another 
possibility  for  validating  models  is  the  higher-order  residual  analysis  of 
Lawrance   and  Lewis    (1986). 
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Joint-maximum    likelihood    estimates   for   k  and    p    (and   perhaps   y) 
obtained   from    the   conditional    density    (2.10)    and   the   formula    for    the    joint 

density  of   X    ,X      ,,•••, X,    is 

n      n-1  1 

f(xn;xn_1;...;x1)    =    fx    |x         <xn;x         )xf  (x         ; x         ) «    . • •    -f      (x    ) , 

n '    n-1  n-1  '    n-2  1 

(2.16) 

where   f„    (x, )    is   a  Gamma(k,1)    density. 
X        1 

Hugus  (1982)  used  (2.16)  to  obtain  joint-maximum  likelihood  estimates 
for  k  and  p.  Three  cases  are  shown  in  Figure  2.1.  Generally,  when  k  is 
greater  than  1,  moment  estimates  of  k  and  p  are  quite  efficient,  which 
agrees  with  results  for  independent  Gamma(k,l)  variates  (Bartlett  and 
Kendall,  1 9 J4 6 )  .  However,  when  k  is  less  than  one,  the  maximum  likelihood 
estimates   become  much   more   efficient    than   the  moment    estimators. 

The  moment  estimators,  p,  k,  p,  given  at  (2.12),  (2.13)  and  (2.14)  serve 
as  good  starting  points  for  numerical  evaluations  to  find  the  maximum 
likelihood  estimates,  u,  k,  p,  of  k  and  p.  Techniques  for  the  numerical 
integration  of    (2.10)    are   given    in  Hugus    (1982). 

3.      MOVING  AVERAGE    AND   HIGHER   ORDER   AUTOREGRESSIVE  STRUCTURES 

The  Beta-Gamma  transformation  can  be  used  to  generate  dependency 
structures  other  than  first-order  aut oregr ess i ve  structures  for  Gamma 
disbributed  time  series.  Several  of  these  structures  are  given  in  this 
section . 

3.1.      The   First-Order   Moving  Average   Process,    BGMA(1). 

The   first-order   moving  average   process,    the   BGMA(1),     is    constructed    in 

essentially    the    same    way    as    the    BGAR(1)    above.      If    {G   }    is   a  sequence   of 

J  J  n 

i.i.d.  Gamma(k,8)  random  variables  and  {B  )  is  an  independent  sequence  of 
i.i.d.  Beta(ka.ka)  random  variables,  where  a  =  1-a,  we  define  the  (backward) 

BGMA(1)  process  { X  }  by 

n 
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Figure  2. 1 

Simulations  of  joint  maximum  likelihood  estimates  for  k  and  p  and  joint  moment 
estimates  for  k  and  p,  for  three  different  sets  of  values  of  the  parameters. 
Ten  replications  for  each  case.  Top  left  figure:  p=0.75,  k=4.0;  top  right 
figure:  p=0.25,  k=4.0;  bottom  figure:  p=0.75,  k=0.75.  The  symbol  (+)  refers 
to  estimates  of  p;  the  symbol  (o)  refers  to  estimates  of  k.  Large  symbols  are 
joint  estimates;  small  symbols  are  marginal  estimates. 
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X   =  G   +  B  G 
n    n    n  n-1 


0  S  a  <  1 


(3-D 


Evidently,  {X  }  is  a  stationary  process  and  X   and  X    are  independent 
n  n      n-r 


>  1.   The  marginal  distribution  of  X  may  be  derived  by  noting  that 


if    \r 

the    right-hand    side    of    (3.1)    is   the   sum   of   a  G(k,8)    random   variable  and  ; 


independent    G(ka,B)    random    variable.       Thus,    X       is    a    G{k(1+a),6)    random 

n 

variable.  This  process  has  the  same  structure  as  the  usual  Gaussian  MA(1) 
process,  except  that  here  the  coefficient,  B  ,  is  a  random  variable  rather 
than    a    constant.       An    immediate    effect    of    this    construction    is   that   the 


observed    and    innovation    processes,     {X    }    and 

n 


G    },    respectively,     have 
n 


different    Gamma    marginal    distributions.       This    is    in    contrast    to    the 

structure  of    the   EARMA   processes    (Lawrance    and    Lewis,    1980),    where    it    was 

deliberately   arranged   that    they   should   have   the   same   distribution.      However, 

as   we  shall    see    shortly,    this    disparity    in   marginal    structure    has    some 

advantages . 

From   the   viewpoint    of    modelling,    it    is   more    useful    to    determine    the 

parameters    of    the    innovation    process    in    terms    of    those    of    the  observed 

process.      For   this   reason,    we  reparameteri ze    (3.1)    slightly.      We    consider 

{G     }     to    be     i.i.d.     Gamma { k / ( 1 +a ) , 6 }     r.v.s.     and     {B    }     to    be    i.i.d. 

Beta{ka/( 1 +a)  ,    ka/(1+a)}    random   variables.      This   yields   an   observed    process 

(X    }    which    is   Gamma(k,B). 
n 

We  may  note  that  if  we  write  p  =  p  (  1  )  =  a  /  (  1  +  a  )  ,  then  G   is 

A  il 

Gamma(kp,6),  the  same  innovation  process  as  for  the  BGAR(1)  process. 

3.2.   Autocorrelation  Function  for  the  BGMA(1)  Process. 

The  autocorrelation  function  for  the  moving  average  process  may  be 
determined  directly.  Thus,  Cov(X  ,X  _  )  =  Cov  ( B  G  ,Gn_  ]  )  =  E(B)Var(G), 
and  so 


P  (r) 
x 


1  +a' 


r|  =  1  , 
rl  >  1  . 


(3.2) 


For  |k|  -  1,  the  attainable  range  of  correlation  is  0  %   p  (1)  i   0.5,  which 
is  the  full  possible  range  of  positive  correlation.   This  is  because,  for  a 
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first-order  moving  average,  |pyM)|  ^  0.5;  see,  e.g.  Hugus  (1982).  The  fact 
that  the  whole  positive  range  is  available  is  important,  because  it  is  in 
contrast  to  the  EMA  models  (Lawrance  and  Lewis,  1977;  1980),  where 
correlation  is  bounded  above  by  0.25.  The  greater  flexibility  in  (3.2)  is  a 
result  of  the  innovation  and  observed  processes  having  different 
distributions.  Since  for  k  =  1,  the  Gamma  distribution  is  an  exponential 
distribution,  the  BGMA(1)  process  is  then  a  broader  first-order  exponential 
moving  average   process    than   the   EMA(1)    process. 

3.3.   Joint  Distributions. 

The  bivariate  Laplace  transform  of  (X   ,  ,X  )  can  be  derived  by  using 

n+ 1   n 

(1.8).   Thus,  again,  using  the  notation  a=1~a,  we  have 


L(u,v)  =  E{exp(-uX   , -vX  }  =  E{exp(-uG   -uB   ,G  -vG  -vB  G   , 
n+1    n  n+1    n+1  n   n   n  n-1 


=  LG(u)E{LG(v+uB)}LBG(v) 


ka  1    „  ikar   „   ika 


B+ul    B+vJ   (B+u+v 


(3-3) 


using    the    Lemma    above.       This    has    exactly    the    same    form    as    the    joint 

transform,    (2.5),    of    (X      , ,X    )    for   the   BGAR(1)    process   with  a  corresponding 

n  +  1      n 

to    p.       This,    too,    corresponds    to  the   behavior    in  Gaussian   processes,    where 

the    joint    distributions    for    the    autor egr essi ve    and    the    moving    average 

processes    have    the    same    form    and    differ    only    in    their    autocorrelation 

functions.      An   immediate   consequence    is   that    the   conditional    distribution  of 

X      ,     given    X    ,    and    X      given    X      , ,    have    exactly    the   same   form   as   for   the 
n  +  1     B  n '  n    &  n+1  '  J 

BGAR(1).      We  note   the   somewhat   unusual    result    for    a    non-Gaussian    process 
that    regression    is    linear    in   both  directions,    even   though   the   process    is   a 

moving  average.      In   fact,    E(X      ,  |X   =x )    =   E(X    |X        =x )    =  ax+ka/8    =  ax+aE(X). 

n  +  1  '    n  n '    n+1 

The    joint    Laplace    transform    of    any    finite    set    of    consecutive 

observations   can   be   obtained   by    the   procedure   that    yielded    (3-3).      Thus,    the 

joint   transform   of    (X    ,X      ,,•••, X  ,)    is   given   by 

n      n-1  n-r+1 


ka    r-1     i      n    ik(1-2a)       <      „     -,ka      r      r  *ka 


L(u, ,u_, • • • ,u    )    =    >- \      x    n      )  {  x    > — —  {      x 

12  r  U+u1  J        i=2    (6  +  u.)  (6+urJ         .=2 

Note    that    this    is    not    the   r-dimensional    transform   for   the   BGAR(1)    process. 

Equality  holds   for    only   r   =   2. 

One     consequence     of     (3.4)      is     that,     since    L(u, u,  •••,._ 

L(u    ,u      ,,•••, u .  )  ,    the   process    is   time-reversible, 
r      r-1  1  r 

The    bivariate    Gamma    distribution  whose   transform   is   given   at    (3.3)    is 
well    known    (Johnson   and   Kotz,    1970b,    p.    219)    and    is   called   by   Ghirtis    (1  I 
the    double    Gamma    distribution.       Since    the    multiv  :  orm    of    this 

bivariate  Gamma   distribution   arises   as    the    individual    sums   of   rn   independent 
Gamma   variates    with   a    common,    independent   Gamma   variate,    it    is   doubtful    that 

triples,    say    X       _,    X      .,    X     ,    in    the    BGAR(1)     process    would    have    tl 

n  +  2    n  + 1    n 

multivariate  distribution.   In  fact,  (3. 4)  shows  that  this  is  not  so  i 

moving  average  process. 

Another  result  that  we  can  immediately  derive  from  the  joint  transform 

is  the  distribution  of  the  sum  of  n  consecutive  observations.   This  has  a 

r 

particularly  simple  form  for  the  BGMA(1)  process.  If  T   =   EX.  then 

n         .     ,     : 

i  =  1 

L   (u)    =    L(u ,u , • • • ,u) ,    which,    from    (3-^),    is   given   by 

k{r-2(r-1 )o} r      „    ik(r-l)a 


Further,    since   0    i    a    S    0.5,    we    can    rewrite    (3.5)     in    terms    of    random 
variables   as 


T      =    G[k(r-2(r-1)a},6]    +    2G{ k (r-1 )a , 6 } , 

where  the  two  Gamma  random  variables  are  independent. 

3.4.   Higher  Order  Moving  Average  Processes. 

Higher  order  moving  average  processes  may  be  constructed  by  extending 
the  BGMA(1)  in  an  obvious  way.   Thus,  the  GBMA(q)  is  given  by 
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q 

X      =    G      +      y      B       .G 
n  n        .    .      n,  1    n-i 

1  =  1 


(3.6) 


where    {G   }    is   a  sequence  of    i.i.d.    Gamma{k/(1+    Ect.),B}    random   variables    and 

1    l 

{B       .}     is    an    independent    sequence    of    i.i.d.    vector   random   variables   with 
n  ,  l 

q  _         q 

B   .,  i  =  1 , 2 ,  •  •  • , q  ,  being  a  {ka. /( 1 +Ea . ) ,  ka./(1+Za.)}  random  variable.   In 
n ,  l  l      l      l      l 

this  case,  {X  }  is  a  stationary  Gamma(k,8)  process  with  X  ,  X    independent 
n  n   n-r 

when  Irl  >  q.   The  autocorrelation  function  is  given  by 


PX(D 


q-1  q 

a   +   )  a  .a .   J  /  (  1  +   /a. 


j-1 


J-1 


0 


r-1 ,2, . . . ,q 


r  >  q, 


which,  again,  is  a  close  analogue  of  the  usual  autocorrelation  function  for 
the  Gaussian  MA(q)  process.  The  major  difference  is  that  all  the 
correlations   are  non-negative. 

4.      THE    MIXED  AUTOREGRESSIVE,    MOVING  AVERAGE    PROCESS,    BGARMA(1,1) 

A  more  complicated  dependence  structure  in  Gamma  distributed  variables 
that    is   the   analog  of    the   usual    linear   ARMA(p ,q )-type   process    is   now   given. 

4.1.      Structure  of    the  Mixed   Process,    BGARMA(1,1). 

We  can  construct  an  ARMA-type  process  with  a  Gamma  marginal  distribution 
by  combining  the  two  first-order  processes  we  have  discussed  above.  For 
convenience,  we  write  each  in  a  slightly  different  form.  The  moving  average 
component    is  given   by 


Y      .    +  B   G    , 

n-1  n  n 


(4.1) 


where    {G    },     (B    }    are    as    given    in   Section    3.1    above,    i.e.,     independent 
n  n 

sequences    of     i.i.d.     Gamma { k / ( 1 +a ) , 8 }     random    variabes    and    i.i.d. 
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Beta{ka/( 1 +a) ,  ka/(1+a)}  random  variables,  respectively.   Notice  that  B   is 

n 

the  coefficient  of  G   in  (4.1),  whereas  it  was  associated  with  G    in 

n  n-1 

(3.1).  Clearly,  this  change  will  make  no  distributional  difference  and,  as 
we  shall  see,  renders  the  parameters  of  the  ARMA  model  more  readily 
interpretable.  The  sequence  {Y  }  is  generated  from  a  BGAR(1)  process  given 
by 


Y   =  A  Y   ,  +  A'G  .  (4.2) 

n    n  n~1     n  n 

The    process    {G   }    is   as    above   and    {A    }    and    {A'}    are   independent   sequences   of 
n  n  n 

i.i.d.    Beta{kp/(1  +a)  ,kp/(1 +a) }    and    i.i.d.    Bet  a { k p/ ( 1 +a  )  , k p / ( 1 +a ) }    random 

variables,    respectively.      If  Y  is  G{ ( k/ ( 1 +a)  ,  8  1  ,    then   so  also   is   Y      and 

the   required    stationary    process    {Y    }    results.       The    product    A'G      is,    of 

n  •  n  n 

course,  the  Gamma{kp/ ( 1 +a) , 8 }  random  variable  used  as  the  innovation  process 
in  Section  2.1  above,  but  is  written  in  this  form  here  to  make  explicit  the 
dependence  on  the  innovation  sequence  G  . 


4.2.   The  Autocorrelation  Function  of  the  BGARMA(1,1)  Process. 

The  autocorrelation  function  of  {X  }  may  be  derived  in  the  usual  way. 

n 

Thus,    Cov(X    ,X         )    =    Cov(Y      , ,Y  , )    +    Cov(Y      ,,B        G        ). 

n      n~r  n-1      n-r-1  n-1      n-r   n~r 

Now,    Cov(Y      , tY  ,)    -    prVar(Y),    and   Cov(Y      , ,B        G        )    = 

n-1      n-r-1  n-1      n-r   n-r 

—  p-  i 

app        Var(G),    so  that   we  obtain 

,     .         p    +   ap      r-1  , .     0> 

PY(r)    =   T~ : P  r   =    1»2,---.  (4.3) 

X  1     +   a 

This  is  the  form  of  the  autocorrelation  function  of  the  ARMA(1,1)  process. 

Note,  toe,  that  the  choice  of  the  slightly  different  structure  for  (4.1)  and 

(4.2)  has  endowed  the  two  parameters  (p,a)  with  physical  significance. 

Choosing  a=0  effectively  sets  B   to  zero,  so  that  the  process  is  now  simply 

n 

p 

the   BGAR(1)    process,    and  we   can    see    from    (4.3)    that    Py(r)    becomes    p    ,    as 

expected.       If    we    choose    p=0,     then    Y      becomes    G      and    {X    }    is   the   BGMA(1 
K  n  n  n 

process,  and  (4.3)  reduces  to  a/(1+a),  as  it  should.  The  reduction  of  the 
mixed  model  to  its  simpler  forms  when  the  parameters  vanish  is  a  consequence 
of    the   structure  we   have   chosen. 
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4.3.      Joint   Distributions    in   the   BGARMA(1,1)    Process 


The  joint   distribution  of    two   consecutive    observations    of    the    process 
} 
Thus , 


(X    }    can    be    derived    in    the    form    of    the    corresponding  Laplace   transform 
n 


L(u,v)    =    E{exp(-uX        -vX    ) 
n+1         n 


E[exp{-u(A    Y        +  A'G    )  '-   uB      ,C      ,    -   vY      ,    -   vB   G   )] 
n   n-1      n  n  n+1    n+1  n-1  n  n 


=    Ln_(u).E{Lv(v+uA)}.E{L_(uA'    +    vB)  }  . 

DO  I  U 


The  first  two  terms  of  this  product  have  already  been  evaluated  and  we  now 
consider  the  third.  By  considering  appropriate  series  expansions,  we  can 
evaluate  the  third  term  in  the  form 


00       00 


I    I 

m=0  n=0 


(-u)m   (-v)n   r(a0+n)r(pQ+m)r(O+m+n)r(e: 


m! 


n! 


r(a0)r(po)r(e+n)r(o+m) 


where  0  =  k/(1+a).   Hence,  we  can  show  that  L(u,v)  is  given  by 


1 


(1+u)(1+v) 


0p 


1 


1  +U+V 


0(a+p) 


^(epsea.e,^), 


(4.5) 


where  F   is  the  Hypergeometric  function,  defined  by 


F  ,       »     y   r(a+n)r(b+n)r(c)    z_ 

2  1 la'D-c'zj    l        r(a)r(b)r(c+n)    n! 
n=0 


The  behavior  and  properties  of  this  function  are  detailed  in  Abramowitz  and 
Stegun  (1964,  Ch.  15).  It  is  easily  verified  that  when  p=0  or  a=0,  the 
appropriate  forms  of  L(u,v)  result  from  (4.5).  The  transform  corresponding 
to  higher-dimensional  distributions  are  more  difficult  to  obtain  in  closed 
form,  although,  in  principle,  series  expansions  of  the  form  of  (4.4)  can  be 
derived . 
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Further,    the   symmetry  of   L(u,v)    in   u  and   v   implies    that   the   conditional 

distributions    of    X  given    X      and    X      given    X       ,     are    identical.        In 

n  +  i         n       n         n + 1 

particular,  we  can  recover  the  conditional  moments  from  (4.5)  and  it  is 
found  that  regression  is  linear  and  the  conditional  variance  is  quadratic  in 
X  . 


4.4.   Higher  Order  Mixed  Processes. 

Higher  order  BGARMA  processes  can  be  derived  by  suitable  extensions  of 
the  BGARMA(1,1).  In  particular,  it  is  straightforward  to  construct  a 
BGARMA(1,q)  process  by  replacing  (4.1)  by  an  MA(q)  form  as  in  (3.6).   Thus, 

q-1 

X   =  I      B   .G   .  +  Y 
n    .  _   n-i  n-i    n-q 
i=0  M 

replaces  (4.1)  and  (4.2)  is  as  before.  The  more  general  problem  of 
extending  to  higher  order  AR  forms  is  more  difficult.  One  way  of  achieving 
it,  however,  is  to  use  mixtures  (random  indexing).  For  details,  the  reader 
is   referred   to  Lewis    (1985). 

5.      DUAL  AND   ITERATED   GAMMA    PROCESSES 

The  first-order  autoregressive  Beta-Gamma  processes  given  in  equation 
(2.1)  has  been  shown  to  be  time-reversible.  This  can  be  a  handicap  in 
modelling  phenomena  such  as  water  run-offs,  which  tend  to  have  'runs  down' 
in  their  sample  paths.  This  is  modelled,  as  noted,  in  a  defective  way  by 
the  GAR(1)  process  of  Gaver  and  Lewis  (1980).  We,  therefore,  look  for  other 
Gamma  processes,  possibly  with  more  than  one  parameter  to  model  dependency 
structure,    which   broaden   the   BGAR(1)    process. 

The  first  process  to  consider  is  the  dual  of  the  BGAR(1).  The  duality 
refers  to  the  fact  that  where  in  (2.1)  we  decrease  the  shape  parameter,  k, 
of    X  by    using    the    Beta-Gamma    transform    and    then    bring    it    up   to   k   by 

adding  an    independent   Gamma   variable,    we  now    increase   k    in  X  by   adding  an 

independent  Gamma  variate  and  then  decrease  the  parameter  to  k  by  using  the 
Beta-Gamma   transform.      Thus,    we  have 
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Xn(k,B)    =   B(k,q){Xn_1 (k,8)    +    Gn(q,B)l.  (5.1) 

However,    it   can   be    shown    that    the    joint    transform    of    X      and    X       ,     is 

n  n 

(1+u)    q(1+v)    ^Cl+u+v)*^        F   {q,q;k+q;uv/(1 +u) (1 +v) } ,    so  that   the   process    is 

time-reversible.      We   thus   have   nothing  new   by   way  of   broadening   the    BGAR(1) 

process. 

Another   approach   to  broadening   the   BGAR(1)    structure    is    to    iterate    the 

process.      Thus,    in    (2.1),    the   left   hand  side    is   a   Gamma(k,8)    random   variable 

and   the   procedure   in    (2.1)    can   be    reapplied.       However,    a    time-reversible 

process    is    again    obtained.      A  better   way   to   iterate   is   to  apply   the   GAR(1) 

procedure    to    (2.1)    and    obtain    a    combination    of    the    BGAR(1)    and    GAR(1) 

processes: 


X    (k,B)    =   Y{B'(kp,kp)X      ,    +    Y    (kp,6)}    +    GI    (kY.B)  (5.2) 

n  n  n-1  n  n 

=    YB'(kp,kp)X      ,    +    YY    (kp.B)    +    GI    (kY,B),  (5.3) 

n  n-1  n  n 

where    0<Y<1,0<p<1,p=Y*1,k>0   and    {GI    (kY.B)}    is   a   sequence 

of    i.i.d.    Gamma   innovation  random   variables   with  Laplace-Stieltjes    transform 

{  (  B+Ys  )/(  8  +  s)  }k,     independent   of    {B'}    and    {Y    (kp,B)l.      The   condition   that    p 

n       n 

and  Y  do  not  both  equal  one  is  necessary  to  obtain  an  ergodic  process.  Note 
that  (5.2)  is  different  from  the  combination  given  in  Lawrance  and  Lewis 
(1982)  and  we  denote  it  by  GBGAR(1). 

Now  in  (5.2),  the  case  Y  =  1  gives  the  BGAR(1)  process,  Y=0  and/or  p=0 
gives  an  i.i.d.  sequence  {X  }  while  p=1  gives  the  GAR(1)  process.  Thus,  we 
should  find  sample  path  behavior  running  from  time-reversibility  to  'runs 
down'  behavior.   Also,  the  process  is  Markovian  and  has  serial  correlation 

p(r)  =  (Yp)  lrl         r  -  0,±1  ,±2,  ••-.  (5.M 


For    the   joint   Laplace-Stieltjes    transform   of   X      and  X      ,  ,    we   have 
J  H  J  n  n-1  ' 


2k 


L       (u,v)  =  E{exp(-uX  -vX  _. )} 


n'  n-1 


E[exp{-uYB^(kp,kp)Xn_1  -  uYYn(kp,6)  -  uGIn(kY,B)  -  vX  _  J] 


B+YuJk(  6      B  lkPr       lkp 


+  u  J  U+uY    8+vJ   [B+v+uY 

Thus,  it  can  be  seen  that  the  process  is  not  time-reversible  unless  Y=1  (the 

BGAR(1)  case).   The  regression  of  X   on  X    =  x  is 

n     n-1 

E(xnlXn~rx)  =  p(1)x  +  {1_P(1^E^X)' 

which  is  linear  in  x,  but  it  is  not  the  same  as  E(X   , IX  =x ) . 

n-1  '    n 

To  separately  identify  the  parameters  Y  and  p  in  this  process,  one  must 
go  beyond  second-order  properties  of  the  process.  This  is  because  the 
parameters  enter  into  the  correlation  (5.4)  as  a  product.  Higher  order 
residual  analyses  (Lawrance  and  Lewis,  1986)  and  maximum  likelihood 
estimation  will    be   considered  elsewhere. 


6.       NEGATIVE    CORRELATION   AND  NON-LINEAR   PROCESSES 

All  of  the  processes  described  above  are  limited  by  their  serial 
correlations  being  non-negative.  There  are  a  number  of  ways  of  extending 
the  processes  to  give  negative  valued  serial  correlations  and  we  discuss  one 
of  them  in  some  detail.  All  methods  involve  non-linear  functions  of,  say, 
X  in    a    first-order   autoregressi ve   process.      This   is   necessary   because  of 

the   non-negativity  and   lack   of   symmetry  of   Gamma   disributed   variates. 

6.1.      Antithetic   Variates. 

Let   X    be   a   continuous    random    variable   with    c.d.f.    F    (x)    and    inverse 

-1  -1 

c.d.f.    F      (a),    0    <    a    <1.       Then    the    random    variable   X*   =    F      { 1 -F    (X)}    is 

A  A  A 

called  the  antithetic  variable  to  X.  For  symmetric  two-sided  random 
variables  centered  at  zero,  X*  =  -X.  For  positive  valued  variables  such  as 
Gammas,    X*    has    the  maximum   attainable    negative    correlation    for    bivariate 
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Gamma  pairs  (Moran,  1967).   In  particular,  if  k  =  1  (Exponential), 

—  x 

X*   =    -£n(1-e      ),    but    if   k*1    the   transformation   is   difficult   to  compute. 

If    X      .in    (2.1)    is    replaced  by   X*    ,    in    (2.1),    then   a   very  non-linear, 
n-1  n-1 

Markovian    first-order    aut or egr ess i ve    process    is    obtained.       Serial 
correlations   beyond   lag  one  are   difficult   to  compute. 


6.2.      Coupl ing. 

Gaver    and    Lewis    (1980)    introduced  a   scheme   in  the   context   of    the   GAR(1) 

processes    for    cross-coupling    two   Gamma    processes    so    that    the  marginal 

processes    will    have   negative   serial    correlations.      It   is   actually  easier    to 

implement    this   scheme   for   the   BGAR(1)    process    than   for    the    GAR(1)  process, 

because  the  random,  Beta  distributed  coefficients  are  continuous.  We  do  not 
pursue   this   here. 

6.2.      Inverse   Processes. 

A  direct  scheme  for  obtaining  negative  correlation  in  a  Gamma  process  is 
now  given.  It  is  a  generalization  of  a  scheme  given  by  Lewis  (1983)  to 
generate  negatively  correlated  bivariate  Gamma  pairs.  Its  utility  lies  in 
the  fact  that  the  sequence  can  be  generated  with  nothing  but  i.i.d.  Gamma 
variates,  no  numerical  inversions  of  inverse  distribution  functions  are 
required . 

Thus,     let    B     (k;q-k),     for    q    >    k,     be    a    sequence    of     independent 

Beta(k;q-k)    variates,     independent    of    G'(q)    and    G"(k+q),    which    are 

v  n    m  n 

independent   sequences   of    independent  Gamma   variates,    n»1  ,2,  •••.       Also,    let 

X^(k)    be   a   Gamma(k)    variate,    where   k   >  0.      The    idea   is   that   we  want   X    (k)    to 
0  n 

be  small  when  X  ,(k)  is  large,  while  retaining  the  Gamma(k)  marginal 
structure  of    the   process.      We   have 

G'(q) 

Xn(k)    =   Bn(k,q-k) (£)+G.(q)    G"(k+q)        q    >    k,    n   =    1,2,....  (6.1) 

n-1  n 

Note  that  the  ratio  is  Beta(q;k)  and,  by  the  Beta-Gamma  transform 
(1.7),    the   product    of    this   ratio  with   the   independent   Gamma(k+q)    variable   is 
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a    Gamma(q)    variate.       The    multiplier   B    (k;q-k)    reduces    the   shape   parameter 
from   q   to  k. 

To  obtain   p(1),    the   correlation   between   X    (k)   and   X      ,  (k),    we  need 

n         n-1 


E{X  (k)X   .  (k)}  =  E{B  (k;q-k)}E 
n    n-1  n   M 


C;(q)Xn.1(k) 


X   ,(k)  +  G'(q) 
n-1       n 


E{C"(k+q) 


_  H  (k*q)  E 
q    B 


G'(q)X   ,(k) 
n    n-1 


Xn-1(k)  +  Gn(q) 


(6.2) 


To  evaluate  the  remaining  expectation  in  (6.2),  we  use  the  fact  that  in  an 

expression  such  as  X   , (k)/{X   , (k)  +  G'(q)},  the  denominator  is  independent 

n-1  n-1  n 

of    the   ratio.      Then   we   have,    after   some  manipulation,    that 


G'(q)Xn_1(k) 


X      ,(k)    +    G'(q) 
n-1  n 


M. 


B(k+q  +  1  ) 


(6.3) 


Combining    (6.2)    and    (6.3).    we  get   the   suprisingly   simple  result   that 


Corr(X    ,X      ,  ) 
n     n-1 


p(1) 


k+q  +  1 


q    >    k 


(6.1) 


This  correlation  is  always  negative  and  if  k  =  1  (the  Exponential  case), 
p(1)  has  a  minimum  attainable  value  of  -1/3  when  q  =  k.  This  is  about 
halfway  to  the  minimum  attainable  correlation  for  bivariate  Exponential 
variables   of   -0.61  . 

The  scheme  can  be  iterated  to  achieve  greater  negative  serial 
correlation.  In  this  and  the  scheme  (6.1),  serial  correlations  of  higher 
order  are  difficult  to  obtain.  However,  since  the  process  is  Markovian,  the 
decay  of  the  absolute  values  of  the  serial  correlations  is  geometrically 
bounded  . 
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7.       AN  ANALYSIS    OF   TIMES   BETWEEN   FAILURES    OF   A    COMPUTER  SYSTEM 

Hugus    (1982)    used   the   Beta-Gamma   process    to  analyze    a    long    sequence    of 

wind    speeds.       This    sequence    is    very    non-stationary,    containing    yearly 

cycles.      The    model    actually    used    is    u(n)G*,    where    y(n)    is    a    log-linear 

function  of   n   and  G*    is   a   unit-mean   BGAR(1)    process. 

n 

A  simpler,  stationary  series  of  times  between  failures  of  a  computer 
system  is  analyzed  here.  Although  this  data  is  known  to  be  generated  by  a 
branching  Poisson  process  (Lewis,  1964),  the  Gamma  model  is  much  simpler  and 
much  more  tractable  than  the  branching  Poisson  process  and  provides  an 
adequate  representation  for  most  purposes.  Modelling  of  these  times  between 
failures  is  important  because,  for  example,  they  represent  times  at  which 
requests   for   service   to   the   computer   are  made. 

In  Figure  7.1,  we  give  a  Gamma  probability  plot  for  the  256  times 
between  failures.  The  fit  appears  adequate,  but  the  goodness  of  fit 
statistics  in  the  table  in  Figure  7.1  must  be  used  with  caution,  since  the 
data  is  serially  correlated  and  two  parameters  have  been  estimated  from  the 
data.  The  parameter  k  (ALPHA  in  the  table)  is  estimated  as  k  =  0.704.  Note 
the  two  outliers  in  the  Gamma  probability  plot.  Since  these  are  serially 
adjacent,    they   probably   represent   a   lapse   in  recording  of   computer   failures. 

In  Figure  7.2,  we  show  a  correlation  analysis  of  the  data.  The  decay  in 
correlation  from  p(1)  =  0.353  to  higher  lags  is  consistent  with  first-order 
autor egr essi ve  correlation  structure.  The  bands  in  the  figure  are 
approximate  confidence  intervals  for  each  p(k)  under  the  assumption  the  true 
correlation  is  zero  for  lag  greater  than  k.  (See  Box  and  Jenkins,  1976,  p. 
35  for  details).  Clearly  the  times  between  failures  are  correlated  and  thus 
a  renewal   model,    say,    for   these  times   between   failures   would   be   inadequate. 

The  partial  autocorrelation  plot  in  Figure  7.2  also  confirms  the  first- 
order   autoregressi ve   nature  of    the   data.      The   last   panel    in   Figure  7.2   shows 

the    autocorrelation    function    for   the   estimated  residuals,    R  =X  — p (  1  ) X      ,; 

n   n      n-1 

there  is  no  significant  correlation  in  this  series. 

Finally,  in  Figure  7.3,  we  give  empirical  density  functions  (kernel 

density  estimates)  for  the  successive  differences,  X  -X   .,  n=2,3i,#*  and 

n      n-1 

successive    ratios    X    /X         ,    n  =  2,3***,    which   were  discussed   in  Section   2.7. 

n      n-1 
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The   differences   show   a   highly   symmetric    and    long-tailed    density    function 

which    is   consistent   with  an    5,-Laplace   distribution.      Note  that   the  median  of 

the   differences    is   estimated   to  be   -9.      Given   the   range  of    the   differences, 

this    is    probably    not    significantly   different   from   the   value  of    zero  which 

would   hold   for   the   time-reversible   Beta-Gamma   process.      The  ratios,    X   /X      ,, 

K  n      n-1 

have    an    estimated   mean    of    10.812    with    estimated    standard    deviation    of 

1  /2 
^2 .  ^32 / ( 255)  =    2.66.       Thus,    there    is    nothing    here    to    suggest    any 

inadequacy    in   the   Beta-Gamma  model    for    characterizing   the   data. 
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Figure  7. 1 

Gamma  probability  plot  for  times 
between  failures  of  a  computer  system 
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Figure  1.2 

Autocorrelation  structure  of  times  between 
failures  of  a  computer  system 
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Figure   7. 3 

Empirical    density   functions    for   ratios    and 
differences    of    successive   times    between   failures 
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